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We consider a model with competing double-exchange (ferromagnetic) and super-exchange (anti- 
ferromagnetic) interactions in the regime where phase separation takes place. The presence of a long 
range Coulomb interaction frustrates a macroscopic phase separation, and favors microscopically 
inhomogeneous configurations. We use the variational Hartree-Fock approach, in conjunction with 
Monte-Carlo simulations to study the geometry of such configurations in a two-dimensional system. 
We find that an array of diamond shaped ferromagnetic droplets is the preferred configuration 
at low electronic densities, while alternating ferromagnetic and anti-ferromagnetic diagonal stripes 
emerge at higher densities. These findings are expected to be relevant for thin films of colossal 
magneto-resistive manganates. 

PACS numbers: 75.47.Gk, 75.30.Kz, 75.10.Lp 



I. INTRODUCTION 

In recent years, doped manganese oxides remained at 
the forefront of theoretical and experimental research 1 . 
The main source of interest in theses systems is the phe- 
nomenon of colossal magnetoresistance (CMR), which 
they exhibit, and that is likely to have important tech- 
nological applications. In the meantime, the underly- 
ing basic physics remains elusive, and probably involves 
the strongly-correlated nature of the doped magnetic ox- 
ides. The CMR in doped manganates is observed for 
intermediate hole-doping levels, typically 0.2 < x < 
0.5, in the temperature region around the transition 
between low-temperature metallic ferromagnetic (FM) 
and high-temperature insulating paramagnetic phases. 
In addition to double-exchange ferromagnetism 2 -, the 
CMR compounds also possess pronounced antiferromag- 
netic (AFM) tendencies, as evident from the AFM spin 
ordering with Neel temperatures of about 100-200K, 
observed 3 ^ at the doping end-points (x = and x = 1). 
This antiferromagnetism is of a superexchange origin 5 . 

The manganates physics involves several degrees of 
freedom of substantially different nature, including lo- 
calized core spins Si of Mn ions, fermionic degrees of 
freedom associated with conduction e g -electrons, lattice 
distortions, etc. In such systems, the presence of compet- 
ing interactions (such as FM and AFM) often gives rise to 
phase separatio n 6 ' 7 ' 8 i 9 , whereby areas of different phases 
are stabilized in a structurally and stoichiometrically ho- 
mogeneous sample. In the case of the manganates, it has 
even been suggested 7 -^ that phase separation into insu- 
lating paramagnetic and metallic FM phases may explain 
the resistivity peak observed near the Curie temperature. 
In the present paper, we focus on the low-temperature 
regime where the presence of phase separation in the 
appropriate manganate systems has been directly ver- 
ified, e.g., by means of scanning tunneling microscopy 



(STM) on thin filmsii. Transport measurements reveal 
metastability and history dependence near the percola- 
tion threshold (x ~ 0.2), confirming phase separation in 
both fil m 12 : 13 and crystallin o 14 : 15 samples. 

Using simple microscopic models^ 6 - it can readily be 
shown— that the hole concentration x indeed controls the 
balance between the FM and AFM tendencies of the sys- 
tem. Once x is tuned away from the optimal CMR doping 
region, the homogeneous FM metallic state no longer cor- 
responds to the energy minimum. Instead, energy can be 
gained by changing the magnetic ordering, carrier den- 
sity, bandstructure, and/or orbital state in part of the 
system, making the sample inhomogeneous 6 ^^ : 16 ' 1 7 . 
The surface tension between different phase a 16 ' 18 then 
competes against the long-range interactions present in 
the system in the form of electrostatic force a 6 ' 16 : 19 ' 20 or 
long-range crystal strain field a 21 : 22 . These require that 
the system remains homogeneous at least on average 
on the appropriate length scale (such as the Debye - 
Hiickel screening length), resulting in a periodic arrange- 
ment of nano- or mesoscopic regions of different phases 6 -^. 
The geometry of the ensuing inhomogeneous (phase sep- 
arated) state is at the focus of our present study. 

Early studies of phase separation in double-exchange - 
superexchange system a 6 ' 16 ' 23 implicitly assumed that the 
effects of the discrete lattice are unimportant, and conse- 
quently treated the problem within the continuum based, 
long- wavelength, approach. Within this framework, the 
optimal phase separation geometry at small values of the 
FM volume (or area) fraction m (also the average mag- 
netization per site) is obviously that of spherical (in two 
dimensions, circular) FM droplets located at the sites 
of a packed hexagonal (triangular) super-lattice. To the 
best of our knowledge, only the three-dimensional case 
was treated in detail, with the implication that in two 
dimensions the situation is similar. When the system 
parameters are varied in such a way that m increases be- 
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yond 1/2, the geometry changes to that of spherical AFM 
droplets in an otherwise FM matrix. The change gener- 
ally occurs via a direct "geometrical phase transition" 23 
without any intervening regime characterized by both 
phases forming infinite connected shapes (such as fila- 
ments and planar slabs in three dimensions or stripes in 
two dimensions ) & 16 i 24 . 

This latter conclusion is important, since such slab or 
stripe arrangements, if realized, would have been charac- 
terized by peculiar and potentially useful properties such 
as history-dependent anisotropy of the ground state re- 
sistivity. However, the continuum treatment, which is at 
the basis of this result, is not valid beyond the region 
of very small values of x <C 1. Indeed, recent studies 
sugges t 17 ' 18 that the boundary between the two phases 
is abrupt on the lattice-spacing scale (i.e., of the type 
commonly associated with Ising spin systems). Such a 
boundary cannot be adequately described in the contin- 
uum limit, and its surface tension depends on its orien- 
tation with respect to the crystalline axes^. This direc- 
tional dependence of boundary energies should in turn 
affect the droplet shape (generally favoring diamond- 
shaped droplets in two dimensions)^, the arrangement of 
droplets in space, and ultimately the way the geometry of 
phase separation evolves with varying m. This is appar- 
ently a generic property of electronic phase separation, 
found also within the frameworks of Falikov-Kimball 25 
and t — J (Ref. H) models. 

In the present paper, we revisit the problem within the 
framework of a single-orbital double-exchange - superex- 
change Hamiltonian (with infinite Hund's coupling Jh), 
augmented by a long-range Coulomb interaction term. 
Using a variational Hartree-Fock approach, we compute 
the energies of various two-dimensional droplet and stripe 
phases corresponding to a FM area fraction m < 1/2, 
and determine the optimal configuration. Our most im- 
portant finding is that while a droplet lattice exists at 
low doping levels, a striped arrangement has a lower en- 
ergy and is therefore stabilized over a broad region of 
the phase diagram. As anticipated from our previous 
results concerning the orientational dependence of the 
boundary energy^, we find that diamond droplets and 
diagonal stripes are the preferred geometries for the FM 
regions of the inhomogeneous states. These conclusions 
gain further support from unrestricted Hartree-Fock cal- 
culations which we have carried out using Monte-Carlo 
simulated annealing on moderate size clusters. The simu- 
lations also demonstrate the existence of inhomogeneous 
states comprised of AFM droplets (or stripes) embedded 
in a FM background (m > 1/2), at higher doping levels. 
While our results pertain to the two-dimensional case, 
it is likely that qualitatively our conclusions would also 
apply to three-dimensional systems. Specifically, we sug- 
gest that a phase separated state with filament or slab 
geometry (rather than a lattice of droplets) is realized for 
a certain range of parameters in three dimensions. 

In addition, we find that the typical droplet size and 
stripe width do not exceed several lattice constants. This 



means that the motion of the charge carriers is strongly 
quantized, rendering droplets midway between metallic 
bulk and magnetic polarons 27 and giving rise to singulari- 
ties in the stripe energy associated with the quantisation 
of the transverse kinetic energy. This important prop- 
erty was not included in the earlier work & 16 ' 23 , which 
assumed sufficiently large length scales for such quan- 
tum effects to be negligible. Our approach, on the other 
hands, allows one to explore the crossover between the 
regime of singly-occupied magnetic polarons, which ap- 
pear for strong Coulomb and AFM couplings, and the 
more conventional phase separation behavior where each 
metallic droplet is populated by many charge carriers. 

The paper is organized as follows: in Sec. HH we intro- 
duce the model and briefly review the physics underlying 
phase separation and magnetic polaron formation in the 
absence of a long-range force. A brief description of the 
calculational methods which were implemented in order 
to include the effects of the long-range Coulomb repulsion 
appears in Sec. Mil while the mass of details is relegated 
to the appendices. Sec. HVl contains a detailed descrip- 
tion of our Hartree-Fock and Monte-Carlo results. We 
conclude with a brief discussion of the results in the con- 
text of current experimental and theoretical work (Sec. 
[Vj ■ While an arrangement of conducting and insulating 
stripes in doped manganate films has not yet been ob- 
served, we suggest that present experimental knowledge 
should allow for a meaningful and successful research ef- 
fort in this direction. 



II. THE MODEL AND ITS PROPERTIES IN 
THE NON-INTERACTING LIMIT 

The starting point for the following calculation is the 
two-dimensional double-exchange Hamiltonian, general- 
ized to include the superexchange coupling and the long- 
range Coulomb interaction, 

n = ~\ £ (4c,a + H.C.) + L ]T S, : ■ Sj 

+ U I, 1 r I ( C » Cm - X ) ( C W - X ) ■ W 

Here t is the nearest-neighbor hopping amplitude and c JQ 
annihilates a conduction electron of spin a =] ,[ at site 
j of a square lattice. Si denotes the core spin made of 
three d-shell electrons (S = 3/2) localized at site i, whose 
AFM superexchange interaction with neighboring core 
spins is given by the second term in H. The third term 
arises from Hund's coupling between the core spins and 
the conduction electrons, where the spin operator for the 
conduction electrons on site i involves the Pauli matrices 
er. It is this term, in conjunction with the fact that the 
hopping preserves the electronic spin, which gives rise to 
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the double-exchange mechanism. This favors a FM spin 
configuration in order to reduce the conduction electrons' 
kinetic energy 2 . The last term includes the Coulomb in- 
teraction among the conduction electrons, whose average 
density is x, and a neutralizing uniform positive back- 
ground, created by the donors. Owing to the long-range 
nature of the Coulomb interaction, the atomic-scale in- 
homogeneities of this background in real systems (cre- 
ated by chemical substitution) are not expected to be 
important from the point of view of our main purpose of 
comparing the energies of various inhomogeneous phases. 
This is because such energies always involve integration 
over volume. 

In using the simplified model, Eq. ([1]), we ne- 
glect some additional physics characteristic of the CMR 
manganates 1 . This includes the presence of two (rather 
than one) conduction electron e g -bands and the electron- 
lattice coupling. The logics behind this simplification is 
summarized, e.g., in Ref. l9j: it is assumed that the mech- 
anism for phase separation (charge ordering) is the com- 
petition between ferro- and antiferromagnetism in the 
presence of a long-range Coulomb repulsion [all contained 
in Eq. |(TJ)]. Once the charge ordering is established, in 
a real system the orbital ordering (and the lattice distor- 
tions) would follow, leading to a quantitative renormal- 
ization of the parameter values. The model, Eq. ([I]), is 
however expected to suffice for a qualitative study of the 
generic features of the phase diagram while its simplicity 
allows to maintain clarity of analysis. Further arguments 
regarding the expected model-independence of our con- 
clusions shall be given in Sec. |V] 

The relatively large value S of the Mn core spins 
means that their fluctuations are small, particularly in 
the T — > limit considered here. In the following we as- 
sume S S> 1 and treat the core spins classically. Conse- 
quently, the effective Hamiltonian governing the physics 
of the conduction electrons is determined by the config- 
uration of the classical spins {S^}. As far as the Hund's 
coupling is concerned the manganates are characterized 
by a moderate bare Jh ^ t. However, they also include 
a strong Hubbard on-site repulsion, Uq S> t, which sig- 
nificantly renormalizes Jh towards the strong coupling 
limit 2 ^. Therefore, while we omit the Hubbard interac- 
tion from our Hamiltonian ([I]), we model its effects by 
taking Jh — ► oo. 

Band theory calculation ;] 9 ' 29 suggest that typical val- 
ues of the hopping amplitude t in the CMR manganates 
lie between 0.3 eV and 0.5 eV. The value of J can 
be roughly estimated from the experimentally observed 
Neel temperatures in the fully doped or undoped (with 
no conduction e g electrons or with no holes) case 3 -^, 
T N ~ 100 - 200K, corresponding to J ~ 5 - lOmeV. 
The long-range Coulomb interaction strength, U , for thin 
films is evaluated as U — e 2 /(aoe). Here, e is the elec- 
tron charge and oo ~ 3.9 A is the lattice spacing. The 
effective dielectric constant e is given by the average of 
dielectric constant e s of the substrate and that of the 
air, e = (e s + l)/2. Among the substances which can 



be used as substrates for manganate films, lanthanum 
aluminate and neodymium gallate hav o 30 ' 31 e s w 23 
and e s w 20, yielding U w 0.31 eV and U « 0.35 
eV respectively. Dielectric properties of the third pos- 
sible substrate, strontium titanate, are strongly depen- 
dent on temperature, wit h 30 ' 32 e s changing from 24000 
(corresponding to U w 0.31 meV) at 4.2K to e s « 277 
(U w 0.027 eV) at 300K. This suggests a possibility of 
experimentally varying the value of U by using different 
substrates and/or changing temperature. 

Theoretical investigations of the double-exchange - su- 
perexchange competition have a history of more than 40 
years. It is by now firmly established^^ that this com- 
petition is resolved not via a second-order phase tran- 
sition from the FM state to a uniform state with a he- 
lical or canted magnetic ordering, but rather via sepa- 
ration of the sample into regions characterized by dif- 
ferent spin arrangements and conduction electron band- 
structures. We will be interested in the case of phase 
separation into FM and AFM regions with abrupt Ising- 
type boundarie o 17 ' 18 between them. This means that the 
resulting configuration of Sj remains collinear, with all 
core spins either parallel or anti-parallel to a selected 
direction. Thus, it is possible to denote a spin state sim- 
ply by Sj/5 = Sj = ±1, and the Hamiltonian of the 
conduction electrons becomes a function of {Si}. The 
large Hund's exchange coupling then forces the conduc- 
tion electrons' spins to polarize in parallel with the core 
spins, resulting^ 3 - in the following distribution of hopping 
amplitudes for a given spin configuration {Si}: 

(rq\\ _ $ —t i,j nearest neighbors and 5j = Sj 
tij {{Si}) - | otherwise 

(2) 

After these simplifications, the Hamiltonian takes the 
form 

« d S i}) = \ E (*« ({Si}) ch + H.c.) + J J2 SiSj 
(*>i) (*>i) 

+ u E fa - *) fa - x ) ■ ( 3 ) 

The model, Eq. ([3]), on the (bi-partite) square lat- 
tice is invariant under the particle-hole transformation 
cj — > (— l) p Ci, and x — > 1 — x, where p takes the values 
0,1 on the two sublattices. As a result we note that in 
the following, x acquire the more general meaning of a 
carrier density, i.e., either the electronic density or the 
hole density relative to the half-filled state. We will now 
briefly review the ground state properties of the Hamil- 
tonian © at U = 0. 

When the carrier density x is finite, the ground state 
of the system at J — > is uniform and FM. With increas- 
ing J beyond a certain critical value J c (x), this uniform 
FM state eventually becomes destabilized, and a non- 
uniform ground state is obtained instead. In this phase 
separated state only part of the system is occupied by 
the FM phase. We will be interested in the case in which 
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the other part is a simple Neel antiferromagnet with zero 
charge-carrier density. A variational studyi2. shows that 
in two dimensions such a phase separated state may be 
realized only for J < J* m 0.036. At higher values of J 
the magnetic ordering in either the electron-rich or the 
electron-poor regions of the sample differs from that of a 
ferromagnet or a Neel antiferromagnet. 
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FIG. 1: The charge carrier density xfm in the FM region of 
a macroscopically phase separated state. 

Thermodynamic equilibrium between macroscopic FM 
and AFM regions means that the thermodynamic poten- 
tials in the two phases are equal, 



where 



and 



Q. FM = I (e - ^)g{e)de + 2 J, 

J — OG 



AFM 



-2J. 



(4) 



(5) 



(6) 



17(e) is the density of conduction electron states in the 
FM region. Since the latter is large, g(e) can be taken 
to be the two-dimensional tight- binding density of states, 
and boundary corrections may be neglected. By solving 
Eqs. ([4][6]) for the Fermi energy, fi, one readily obtains 
the carrier density in the FM region, xfm, via 



xfm 



g(e)de, 



(J) 



with the result depicted in Fig. [TJ The system remains 
in a uniform FM state as long as xfm(J) < x. The 
fraction of the system area (or volume) , occupied by the 
FM phase, is given by to = x/xfm- The critical value, 
J c , for the onset of phase separation is then determined 
by the condition Xfm(Jc) = x - 



In addition to the macroscopic phase separation as de- 
scribed above, the double exchange-superexchange com- 
petition can also be resolved via an altogether different 
scenario (formation of magnetic polarons). When a sin- 
gle electron (or hole) is lodged into an antiferromagneti- 
cally ordered double-exchange magnet with zero charge- 
carriers, (x = Q), a free (self-trapped) magnetic polaron, 
or ferron£ i 27 i 34 i 35 i 36 i 37 ' 38 , is formed around it. It is es- 
sentially a microscopic FM region, containing one charge 
carrier, in an otherwise AFM system. Since the propaga- 
tion of charge is unimpeded in the FM region it acts as 
a potential well for the sole carrier, which occupies the 
lowest bound state inside the well. The polaron binding 
energy, E mp , (with respect to the state where the AFM 
order is unperturbed) can be be easily estimated^ 7 -. We 
consider the case of a diamond-shaped polaron , with 
L + l sites along each side (see Fig. [2 upper left). For 
L> 1 we find 



E mp (L) 



-It- 



L + l 



8L 2 J, 



(8) 



where the first two terms are the ground-state energy of 
the charge carrier and the last one represents the superex- 
change contribution. Expression ijHJ should be minimized 
with respect to L, resulting in 



Emp — 



-2< + 4ttV Jt 



(9) 



Here, the coefficient of the second term depends on the 
geometry of the FM micro-region (e.g., for a round po- 
laron one would have obtained 12.06 instead of iw). 

The above expressions are valid in the J/t <C 1 regime 
(yielding L ~ (i/J) 1 / 4 S> 1), where it is easy to verify an 
important statement which is expected to hold for all J. 
Namely, if in the absence of a Coulomb interaction U, a 
second carrier is added to the system, it is energetically 
favorable for the two charge carriers to occupy the two 
lowest bound states in a shared FM micro-region rather 
than to form two independent polarons. This conclu- 
sion is verified by calculating the binding energy of the 
(diamond-shaped) doubly-occupied polaron 



mp 



-At + 4ttv/7J</2, 



(10) 



which clearly satisfies Emp < 2E mp . This trend continues 
when further charge is added, and at n 1 the binding 
energy (per carrier) of the n-carrier polaron decreases 
toward the limiting value E ps , 



l E {n) 
mp 



E„ 



(11) 



which is the energy gain per carrier associated with the 
macroscopic phase separated state. The latter can be 
evaluated as 



En 



Efm — Eafm 



XFM 



(12) 
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where Efm{h) = JT^ sg(e)de + 2J and Eafm = — 2J 
are the energies per site of the FM and AFM phases. 
Using Eqs. (j4|) and {2$ to evaluate [i and xfm we find, 
in the limit J<t, 



(13) 



The inequality, E ps < E mp , implies that for any finite 
carrier density, at U = 0, the double-exchange - superex- 
change competition is resolved via macroscopic phase 
separation. 

Notwithstanding the preceding discussion, its conclu- 
sion may change if a realistically strong Coulomb inter- 
action U is included, favoring a large spatial separation 
between the charge carriers. Indeed, as we demonstrate 
in the following, a polaronic state arises in the regime of 
large U and small carrier density. It is the extreme limit 
of a broad range of inhomogeneous states which originate 
from the frustration of macroscopic phase separation by 
long-range forces. The study of this intermediate region 
of parameters lies at the focus of the remaining part of 
the paper. Since the typical size of the resulting FM 
regions is rather small one needs to take into account 
the effects of quantization of the charge carrier motion. 
At the same time, some of the results obtained for the 
macroscopic phase-separated system, such as the direc- 
tional dependence of the boundary energy^, still offer 
important guidance to the understanding of the inhomo- 
geneous configurations. Next, we outline the methods 
used to treat this intermediate regime which is charac- 
terized by a combination of both traditional phase sepa- 
ration and magnetic polaron (quantized) physics. 



III. METHODS 
A. Variational Hartree-Fock Approach 

Given the Hamiltonian , our task is to find the con- 
figuration of core spins in the ground state. However, 
there is a vast multitude of possible spin configurations 
amongst which the ground state is to be sought, making 
it impossible to explore all of them. Nevertheless, previ- 
ous studies of similar or related systems^ suggest several 
families of highly symmetrical configurations as reason- 
able ground state candidates. The two main types of 
spin configurations studied in this work are FM droplets 
in an AFM background, and alternating FM and AFM 
stripes, as illustrated in Fig. [2] A uniform FM phase, in 
which the double-exchange mechanism completely over- 
comes the superexchange, is also considered. 

Calculating the energy of the conduction electrons in 
a given configuration of core spins is a difficult problem. 
Here we suffice with the Hartree-Fock (HF) approxima- 
tion, which gives an upper-bound to their ground state 
energy. Since we are dealing with periodic spin config- 
urations, the HF equations for the whole system can be 
rewritten as an effective eigenvalue problem within a sin- 
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FIG. 2: The geometry of the inhomogeneous configurations 
considered in the variational Hartree-Fock calculation. The 
dashed lines indicate the super-lattice unit cell. 



gle unit cell. The superexchange contribution to a con- 
figuration's energy is simply calculated by counting the 
number of FM and AFM bonds in a unit cell. 

Based on previous analytical results^ and numerical 
investigations^ the considered droplets are either dia- 
mond or square shaped, and are chosen to form either a 
triangular or square super-lattice (see Fig. [2]). Several 
droplet phases are possible by combining different droplet 
shapes and super-lattice types. In addition, one has vari- 
ational freedom to specify L, the size of the FM droplets, 
and n, the number of conduction electrons in each one 
of the droplets. The distance between the droplets, a, is 
then uniquely determined by the type of super-lattice, by 
n, and by the average density of conduction electrons, x. 

The energy of each droplet phase is found by minimiz- 
ing its energy density, 



E, 



droplet 



Ehf ^droplet | j-i /-i a \ 

X h rj J, droplet ; (14) 



with respect to the variational parameters, L, and n. 
Here EnF,dropiet is the HF energy of the conduction elec- 
trons inside a unit cell containing a single droplet, and 
Ej^dropiet is the AFM energy per site. The details of the 
HF calculation appear in Appendix El The main source 
of complication is the necessity to take into account the 
Hartree interaction between electrons belonging to dif- 
ferent droplets in the infinite super-lattice. This is done 
by employing Ewald's summation method (see Appendix 
[Bl) . The AFM coupling energy per lattice site is 



E 



diamond droplets 



-2 J ( 1 - 5L ~~ 4 8£+8 ) square droplets 



J, droplet 



(15) 

where A ~ n/x is the area of a unit cell. 

Two types of stripe phases were considered: diagonal 
and bond-aligned. Additional variational freedom comes 
from the need to specify W, the stripe width, and xfm, 
the (average) density of conduction electrons within the 
FM stripe. Just as for the droplet phase, the stripe phase 
energy is found by minimizing its energy density 



E 



Ehf. 



stripe 



stripe 



X 



E 



J, stripe 



(16) 
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with respect to W and xfm- Here A and Ehf, stripe are, 
respectively, the conduction electrons' number and en- 
ergy per unit cell, and Ej tStr i P e is the AFM energy per 
site. A unit cell in diagonal stripes is only one lattice 
spacing long in the direction along the stripe, and two 
spacings long in bond-aligned stripes (see Fig. [21); its 
width equals the stripe periodicity. Therefore, 



XfmW diagonal stripes 
2xfmW bond-aligned stripes 



(17) 



A, together with x, uniquely determine the distance be- 
tween stripes D. 

In a similar manner to the case of the droplet phase, 
when calculating the HF energy one needs to take into 
account the Hartree interaction between the infinite num- 
ber of unit cells in the systems. Moreover, the extended 
nature of the states along the stripes means that it is 
necessary to consider also the Fock exchange between 
different unit cells on the same stripe. A detailed ac- 
count of the way this is done is presented in Appendix 
lAl The AFM spin coupling energy per unit area for both 
diagonal and bond- aligned stripe phases is 



E 



J, stripe 



-23 



1 -2m 2 



1 

W 



(18) 



where m = x/xfm is the fraction of FM regions in the 
system. 

By comparing the energies of all the above mentioned 
phases, a phase diagram is constructed, depicting the 
nature of the ground state as a function of the external 
parameters, x, 3 ft and U/t. 



and incrementally flipping the spins on the other sub- 
lattice. An additional improvement comes from a new 
algorithm used to decide which spin to flip. At first, a 
spin is chosen randomly. It is flipped if the resulting state 
satisfies the Metropolis criterion. If the spin is located 
near a FM-AFM boundary then its neighbors are added 
to a queue of spins to be tested for flipping. After all 
the spins in the queue have been tested for a flip, a new 
spin is chosen randomly. Requiring that a spin be added 
to the queue no more than once, prevents the simulation 
from repeating itself, thus maintaining ergodicity. 

In our calculations, the system contained 24 x 24 sites 
arranged periodically on a torus. A linear annealing 
schedule was employed over 50 — 100 MC sweeps, and 
an identical number of sweeps at the lowest temperature 
allowed the system to thermalize into the ground state. 
The temperatures started from above 1.53 at the begin- 
ning of the annealing schedule to below 0.5 J at its end. 

Even though this method minimizes the system's en- 
ergy with respect to an unrestricted configuration space, 
it has a number of disadvantages when compared to the 
variational HF method, applied to only a number of spe- 
cial configurations. First, its periodicity is fixed; in our 
case it is 24 sites along each axis. In addition, the pres- 
ence of long range interactions causes the simulations to 
converge very slowly. Nevertheless, it provides an impor- 
tant reference point with which the the variational HF 
results may be contrasted, especially in order to confirm 
that the variational manifold contains the most relevant 
configurations. 



IV. RESULTS 



B. Monte-Carlo Simulated Annealing 

We have supplemented the calculation of the HF en- 
ergy of various variational configurations by an unre- 
stricted HF calculation using Monte-Carlo simulated an- 
nealing. In this method, the energy of a finite-sized sys- 
tem is minimized with respect to the full configuration 
space of core spins, rather than a special subset of spin 
textures. In each Monte-Carlo step the energy of a given 
configuration of classical core spins is evaluated using the 
HF approximation. A spin configuration is accepted as 
the system's new state if the change in energy from the 
current state satisfies the Metropolis condition. The tem- 
perature is slowly decreased until the system reaches a 
stable, low energy configuration. If the temperature is 
decreased slowly enough, the final state is the HF ap- 
proximation of the ground state. 

The underlying assumption of the present study is that 
the system indeed separates into FM and AFM regions 
with an abrupt boundary between them. Therefore, the 
MC simulation needs to explore only such configurations, 
improving the convergence time. This can be achieved by 
setting all the spins on one sub-lattice to the "up" state, 



In the present section, we present and analyze our nu- 
merical results. The coupling constants 3 and U are 
measured in units of t, by setting the hopping amplitude 
t = 1 . The HF energies of all the considered phases were 
calculated in the parameter range x < 0.1, 3 < 0.03, 
and for three values of U, namely U = 0.025, 0.075, 0.25. 
We chose to concentrate on this region in the x — 3 
plane for two reasons. As already mentioned, a previ- 
ous estimate^ sets 3* = 0.036 as the upper limit for 
the realization of a FM - Neel AFM (as opposed to other 
types of magnetic ordering) phase-separated state in two- 
dimensions. Secondly, our calculations indicate that the 
line x/xfm = 1/2, crosses 3 — 0.03 at x* = 0.1, see Fig. 
[51 The region below this line in the x — 3 plane corre- 
sponds to configurations in which the FM phase occupies 
more than half the system area. While the stripes phases, 
which we consider, continue to be relevant in this region, 
we expect (and confirm in our MC simulations) that the 
phases of FM droplets ought to be replaced by configu- 
rations of AFM droplets embedded in a FM background. 
The latter turn out to be more involved computationally 
and were left out of the present study. We also wish to 
note that the above values of 3* and x* , are sensitive 
to the details of the considered model. Therefore, while 
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FIG. 3: The HF phase diagram: DD - A triangular lattice 
of diamond-shaped droplets, SD - a square lattice of square 
droplets, DS - diagonal stripes, FM - a uniform ferromagnet. 
The black lines correspond to x/xfm = 1/2. 

experimentally, percolation of the metallic phase at low 
temperatures is observed in manganates with x > 0.18, 
we expect our qualitative conclusions to apply to more 
complicated models of manganates, as long as phase sep- 
aration into FM and AFM phases is possible. 

We begin our review of the results by discussing the 
phase diagram and presenting general arguments for its 
structure. We then move on to consider the details of the 
most dominant phases. 

A. The Phase Diagram 

The main result of our calculation is the phase dia- 
gram, Fig. [31 derived from the variational HF approach 
and depicting the system's ground state configuration as 
a function of the parameters, x, J and U . It demon- 
strates that diamond droplets in a triangular formation 
is the preferred phase at low densities, while diagonal 
stripes are prevalent at higher values of x. This stability 
of a striped phase is the most important result of our 
calculation. The striped arrangement is expected to pos- 
sess unusual and potentially useful properties (see Sec. 
[Vl below, where we also mention possible directions of 
experimental search for the stripe phase in CMR man- 
ganates). When the Coulomb interaction strength U is 
increased, the transition between droplets and stripes oc- 
curs at higher x. As discussed above, our variational ap- 
proach becomes insufficient below the line x/xfm = 1/2, 
as we do not allow for a phase of AFM droplets embed- 
ded in a FM background. Such a phase is expected to 
appear near the transition to the uniform FM state. This 
conclusion is supported by the unconstrained HF results 
presented below. We are unable, though, to map in de- 
tail the boundary between the stripe and droplet phases 
in this parameter regime. 

The general features of the phase diagram can be ex- 
plained by simple energy considerations. The preference 
of diagonal stripes and diamond droplets is a direct result 
of the lower energy of diagonal boundaries, as previously 
established by two of the authors^. The appearance of 
a triangular droplet lattice at low densities is akin to the 
physics giving rise to the Wigner crystal in a dilute gas of 



electrons. Next, we elaborate on the reasons and nature 
of the transition between the droplet and stripe phases. 

To this end, let us examine how the energy difference 
between the two phases evolves with x. As x increases, 
the distance between droplets or stripes diminishes, but 
our HF results indicate that the size of these FM regions 
and the electron density within them, xfm, vary slowly 
in the vicinity of the phase transition. The combined 
difference between the kinetic and magnetic energies per 
electron of the two phases, Ae, depends on xfm and 
on the size and shape of the FM regions, but not on 
the distance between them. We therefore conclude that 
at the qualitative level, changes of Ae with x cannot 
be the driving force behind the transition. Instead, we 
concentrate on the doping dependence of the difference 
in the Coulomb energy per electron between droplets and 
stripes, A<f>. 

The Coulomb energy contains contributions coming 
from the interaction between charges within a single 
super-lattice unit cell and between different cells. The 
neutrality of each unit cell (due to the positive back- 
ground) implies that the dominant contribution to the 
Coulomb energy per electron originates from the intra- 
cell component. Simple dimensional analysis allows us 
to obtain an estimate for its behavior. The amount of 
positive background charge within a droplet unit cell 
is xa 2 , a being the inter-droplet spacing. Thus, the 
Coulomb potential due to the positive background is 
^droplet ~ —Uxa. The interaction between electrons 
within a droplet generates 4 > dropiet ~ UxfmL, where L 
is the droplet size. Since m = x/xfm ~ L 2 /a 2 , we have 

4>droplet = 4 droplet + ^droplet ~ Ux FM L (l - y/m) . (19) 

The Coulomb energy in the stripe phase takes a dif- 
ferent form. The amount of charge per unit length is 
xfmW = xD where W and D are the stripe width 
and the distance between stripes, correspondingly. The 
background potential is then 4>g tripe ~ UxfmW In D 
and the potential due to electrons in the same stripe is 
4>7tri P e ^ -UxfmW In W. Together they give 

4>st n pe = 4>7t n pe + 4>ttripe ~ ~ U XF M W In TU, (20) 

where in this case m = W/D. Consequently, the differ- 
ence in Coulomb energy per electron between the droplet 
and stripe phases has the form 

A0 = 4* droplet 4*stripe 

pa Ux FM [K d L (1 - ^/m) + K s W\nm] , (21) 

where Kd and K s are numerical constants character- 
izing the geometry of droplets and stripes, respectively. 
It is implicitly assumed in Eq. (|2ip that in the transi- 
tion region between the phases Xfm is the same for both 
configurations (the HF calculation shows that this is cor- 
rect up to 20%). The transition itself takes place at m*, 
satisfying 

A0(m*) + Ae = O, (22) 
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FIG. 4: Dimensional analysis form of A(j> as a function of 
m, demonstrating the difference between the two possible 
branches. 



the orientational dependence of the boundary energy and 
the quantization of the carrier motion.) 

Fig. [5] shows HF results for Acf> and (— Ae) as a func- 
tion of x for various values of J. Transitions between 
droplet and stripe phases occur when A<f> = — Ae. The 
division into discontinuous segments is due to changes in 
the properties of the stripes or droplets respectively (the 
optimal values of W, L, and xfm, see below). However, 
the transitions generally do not occur at these points of 
discontinuity, leading to our previous assertions concern- 
ing the constancy of Ae and the dominant role of the 
Coulomb interaction in the vicinity of the transition. The 
first two transitions at J = 0.014 and J — 0.0165 are near 
a maximum in A</>, demonstrating the 2K S W < KdL 
branch behavior. Whereas these transitions occur on one 
continuous segment, the third transition occurs on a dif- 
ferent segment, where only one solution exists. Other 
one solution transitions are shown for J = 0.0225 and 
J = 0.0275. 
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FIG. 5: A(p and — Ae as a function of x as obtained from 
the HF calculation for U = 0.075 and various values of J. 
Below each plot are colored bands showing the ground state 
configuration at the corresponding x values: green (light)- 
droplets, blue (dark) - stripes. The transitions occur when 
Acj> = -Ae. 



B. The Droplet Phase 

In general, a triangular lattice of diamond shaped 
droplets proved to be energetically more favorable than 
the other types of droplet phases. As noted before this 
is a consequence of the directional dependence of the 
FM-AFM boundary energy and the minimization of the 
inter-droplet Coulomb energy. Fig. [S] shows the optimal 
droplet size, L, and the number of conduction electrons 
per droplet, n, as deduced from the variational HF cal- 
culation. Increasing the strength of the Coulomb repul- 
sion has the obvious effect of decreasing the droplet size. 
Specifically, for the case of U = 0.25 the variational study 
yields very small (L = 1,2) singly occupied droplets in 
the regime of low x and intermediate to large J. Com- 
paring their energy to the other types of inhomogeneous 
states reveals that these magnetic polarons are in fact the 
lowest energy configuration in this region of parameters, 
see the HF phase diagram, Fig. [3] 




where we have used the constancy of As near the tran- 
sition. 

If 2K S W > KdL then A<^> is a monotonously increas- 
ing function of m (in the physical range < m < 1), see 
Fig. |U In this case at most a single solution, to*, exists 
to condition (|22"|) , implying that the droplet phase is pre- 
ferred when to < m*, while stripes occur for to > to*; 
the area of the droplet phase increases with U. On the 
other hand, if 2K S W < KdL, A<f> acquires a maximum 
and two solutions, to* and m*>, may appear. Under such 
conditions a reentrant behavior follows, i.e., droplets are 
preferred when to < to* or to > to?,, and stripes are 
realized in the region to* < to < to?., which grows with 
increasing U. We note that in any case the existence of a 
solution to condition (f2"2"]l , crucially depends on the value 
of Ae. It is the latter which reflects the features taken 
into account for the first time in the present work (viz., 
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FIG. 6: HF results for a triangular lattice of diamond 
droplets: optimal size L (top) and number of electrons per 
droplet n (bottom). 
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C. The Stripe Phase 

Fig[7| shows the optimal stripe width W and conduc- 
tion electron density Xfm f° r diagonal stripes. The lat- 
ter are more favorable than their bond-aligned counter- 
parts due to the orientation dependence of the boundary 
energy. One striking feature in these HF results is the 
existence of abrupt transitions in the stripe width. A 
small increase in x may lead to a discontinuous change 
in W . On the other hand, increasing J typically leads 
to changes in W which are less steep. The electron den- 
sity within the FM stripes, xfm, varies, in general, very 
slowly with x, and increases with J. 



U = 0.025 



U = 0.075 



C/ = 0.25 
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FIG. 7: HF results for diagonal stripes: optimal W (top) and 
xfm (bottom). The number of partially filled bands within a 
stripe changes across each black contour on the top panel. 

We use the condition of thermodynamic equilibrium, 
Eq. ([4|, between a diagonal FM stripe and its AFM en- 
vironment to explain these features. The kinetic energy 
contribution to the stripe's thermodynamic potential is 
determined by its non-interacting spectrum consisting of 
W bands (corresponding to the quantization of trans- 
verse electron motion within the stripe) 



£b{k) = -ifceos ( - 



where % is the bandwidth of band b = 1 • • • W 

&7T 



tt, = 2icos 



W + l 

The resulting density of states in band b is then 

/ \ 2 1 
9b(s) 



■nth 



1 - (e/t b y 



(23) 



(24) 



(25) 



which together with the chemical potential fi determines 
the number of electrons rib per unit length in the band. 
Since we are interested in relatively low doping levels we 
consider the lower W/2 bands for which 



n b 



n 2 

gb{e)de = — cos" 
-n * 



and the non-interacting electronic contribution to the to- 
tal energy is 



E b = 



egb(s)de 



2t 



Using these terms, the thermodynamic potential in the 
FM stripes is 

n FM (w) = i J2( E * - + 2J f 1 - |f) . (28) 



where the second term is the magnetic energy, taking into 
account the structure of the boundaries. £Iafm remains 
the same as for an infinite AFM region, Eq. Fig. [5] 
shows xfm — J2b n b/W evaluated at the Fermi energy /j, 
which solves Qfm(W) = £Iafm as a function of J and 
W. 



xfm 




(26) 
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FIG. 8: xfm in diagonal stripes of width W as imposed by 
thermodynamic equilibrium. 

Points of non-analyticity occur whenever the chem- 
ical potential increases beyond the bottom of a band, 
fi = —tb, so that the carriers begin to fill this additional 
band. These non-analyticities result in a corrugated 
landscape for xfm(W, J), shown in Fig. [SI whereby sev- 
eral values of W may correspond to the same xfm- Thus, 
a small increase in x may drive an abrupt change in W 
but leave Xfm constant. This transition is accompanied 
by a change in the number of partially filled bands within 
the stripe, as depicted by the black contours in Fig. [7] 



D. Simulation Results 

Although our Monte-Carlo results are not sufficient for 
constructing the phase diagram, they yield convincing 
evidence that the phases included in the variational HF 
calculation are indeed the appropriate variational phases 
to consider. Some examples of ground states obtained by 
MC simulated annealing are given in Fig. [Sj Note that 
the moderate cluster size used in the simulation induces 
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finite-size effects apparent, for example, in the imper- 
fections of the stripe configurations. These results, in 
addition to results from other simulations done at other 
parameter values, agree with the general structure of the 
phase diagram in Fig. [3] Moreover, the unrestricted na- 
ture of the MC method yields also configurations with 
AFM droplets in a FM background (bottom right of Fig. 
[5]) . As mentioned before, such states were not considered 
in the variational HF approach because of the relative 
difficulty in calculating their HF energy. Nevertheless, 
there are reasons to believe, as is confirmed by the simu- 
lations, that such a phase indeed exists around the tran- 
sition line between the striped and uniform FM phases, 
where m — x/xfm > 1/2. 




FIG. 9: (color online) Electron density in ground states 
reached using simulated annealing. Top left: x = 0.0278, 
J = 0.015, U = 0.075. Top right: x = 0.0486, J = 0.02, 
U = 0.05. Bottom left: x = 0.0556, J = 0.015, U = 0.05. 
Bottom right: x = 0.0833, J = 0.01, U = 0.05. Dark lines 
outline FM-AFM boundaries. 



V. DISCUSSION 

The main finding of this paper concerns the geometry 
of the low temperature phase-separated state in a two- 
dimensional double-exchange magnet. We did not in- 
voke any lattice or orbital degrees of freedom but instead 
concentrated on the effects of the ubiquitous long range 
Coulomb interaction. We verified that when the relative 
area occupied by the FM phase, x/xfm, is sufficiently 
large, a striped arrangement (rather than a droplet super- 
lattice) is stabilized. Our results also confirm the expec- 
tation, based on a previous analysis of the directional 
dependence of the FM-AFM boundary energy-i 8 , that 
diamond-shaped droplets and diagonal stripes are pre- 
ferred over their square and bond-aligned counterparts. 
Indications to this effect are also present in other numer- 
ical studies of double-exchange models^. 

The stability of the stripe phase should not come as a 



surprise. In fact, even in the earlier studies, which consid- 
ered the continuum limitSiS, it was noted that the ener- 
gies of the stripe and droplet configurations can be very 
close, although no parameter window was found where 
stripes would correspond to the lowest energy configura- 
tion. As a result, it is conceivable, as indeed was shown 
in Ref. 40i, that a stripe phase may be stabilized, even in 
this limit, once physics due to some sort of additional de- 
grees of freedom is taken into account. Stripes also occur 
naturally in other models, such as t — J or Hubbard^, 
which involve a competition between the AFM nature 
of a parent undoped state and the kinetic energy of the 
doped charge carriers. Long-range AFM interaction was 
found to favor stripes in the FM Ising model^ 2 -. Re- 
garding the case of a pure double-exchange system with 
Coulomb repulsion considered here, it has already been 
argue d 18 ! 20 that the correct treatment of the boundaries 
between the FM and AFM regions is likely to tilt the 
balance in favor of a striped arrangement. 

In the present work, we considered the experimentaly 
relevant case of nanometer-size FM inclusions (compris- 
ing only a few lattice periods). Beside mapping the evo- 
lution of the geometry of the inhomogeneous system, we 
addressed the long-standing question regarding the sta- 
bility of free magnetic polarons. As expected, we find 
polaronic behavior in the region of small carrier concen- 
tration x <C 1 and strong Coulomb interaction. Away 
from this regime, individual magnetic polarons coalesce 
into larger FM areas. We were able to span the entire 
intermediate regime between the coventional phase sepa- 
ration (where the quantized character of the carrier mo- 
tion becomes unimportant) and an array of free magnetic 
polarons (for which the notion of thermodynamic equilib- 
rium between FM and AFM phases becomes irrelevant). 
We emphasize that the two main physical ingredients un- 
derlying our findings, namely, the quantized electronic 
motion in small FM regions and the directional depen- 
dence of the boundary energy, can be viewed as largely 
model- independent. Therefore, our present conclusions 
can be expected to stand for any double-exchange model 
with a long range interaction, including the case when 
the latter originates from crystal strain field a 21 ' 22 . 

The bulk of our study was carried out using a varia- 
tional HF approximation for the energy of various droplet 
and stripe phases. It was supplemented by unconstrained 
HF calculations on moderate size clusters, implemented 
via Monte-Carlo simulated annealing. The HF approxi- 
mation is expected to gain accuracy whenever the ratio of 
electrostatic energy to kinetic energy is small. Through- 
out the range of parameters studied by us, this ratio never 
exceeds 0.15. Moreover, since we deal with the case of 
fully polarized electronic spins, the spatial part of the 
many-body wave-function is antisymmetric. This fact 
reduces correlation corrections to the HF result which 
stem from the tendency of any pair of electrons, owing 
to their mutual repulsion, to be more distant from each 
other than the HF wave-function would indicate. 

We close with a brief discussion of the experimental sit- 
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uation. To the best of our knowledge, a conclusive experi- 
mental observation of metallic stripes in phase-separated 
films of CMR materials is yet to be made. We note, 
however, that stripe-like charge ordering on the atomic 
scale (charge density wave) was observed in a variety of 
manganates. This includes films with different doping 
levels^ 3 -, as well as cerami&M and single crystal samples 
in the insulating state above x = 0.5. In addition, it was 
suggested^ 6 - that the phase separated state in a three- 
dimensional system may acquire a filament structure. 

Nevertheless, we argue that it would be desirable to 
synthesize manganate films whose phase-separated state 
clearly exhibits metallic stripes. In addition to illustrat- 
ing our theoretical picture, such systems are expected 
to display unique and potentially useful properties, some 
of which were not previously observed. One of these is 
an anisotropic conductance, whereby the stripes direc- 
tion determines a low-resistivity axis, which ought to be 
amenable to reorientation by, e.g., applying a voltage. In 
general, one expects to find history-dependent resistance 
and memory effects akin to, and probably more pro- 
nounced than those observed earlier in phase separated 
film o 12 ! 13 , for which not evidence for stripes was reported. 
When the sample composition gets close to the one corre- 
sponding to a stable striped arrangement, weak perturba- 
tions such as external electric or magnetic fields may be 
sufficient to change the geometry of the FM regions from 
droplets to stripes, with a drastic change in transport 
properties in the form of colossal electroresistance due 
to dielectrophoresis^i and large low-temperature magne- 
toresistance. 

Which manganate system could potentially exhibit a 
metallic stripe order? In general, in order to look for such 
a state one is interested to explore the parameter space by 
changing the average carrier concentration x, the metal- 
lic area fraction m, and the strength U of the Coulomb 
interaction^. While x is determined by the dopant con- 
centration, the experimentally measurable quantity m 
depends, in our model, on the ratio of the AFM cou- 
pling J to the hopping t. The latter may be affected by, 
e.g., the choice of the rare earth ion. An example of a 
system which apparently allows control over the value of 
m is (Lai_ y Pr y )i_ x Ca x Mn03. Three-dimensional crys- 
tals of this compound with x between 0.25 and 0.5, are 
metallic for y = 0, with no signatures of phase sepa- 
ration at low temperatures 3 . At y = 1, the system is 
phase-separated^, and exhibits robust insulating behav- 
ior, presumably corresponding to well-separated metallic 
droplets in an insulating matrix. The properties of the 
phase-separated state change as one decreases the value 
of y, and at y = 0.7 it is possible to observe conduction 
paths formation and switching as a result of an applied 
current 15 . Similar behavior is also found in thin films of 
the same compound, which at least for sufficiently large 
values of y are phase-separated^, as reflected in their 
peculiar dielectric and transport propertie a 51 ' 52 . 

These findings prompt us to suggest looking for signa- 
tures of stripes in (Lai_ y Pr y )i_ x Ca x Mn03 films by sys- 



tematically varying y, and with it, as indicated above, 
the relative area m of the metallic phase. We expect 
stripes to appear around the point where the areas of 
metallic an insulating phases are equal to each other. 
Besides (Lai_ y Pr v )i_ x Ca x Mn03, there are other hole- 
doped manganate systems which may exhibit a stripe 
geometry of phase separation, see Ref. l53l In addition, 
we expect our results to be relevant for some electron- 
doped manganates 5 ^, as well as possibly for Eu-based 
magnetic semiconductors 6 . 
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APPENDIX A: HARTREE-FOCK EQUATIONS 
FOR A PERIODIC CONFIGURATION 

The HF equations for n interacting, spin polarized, 
electrons, may be written in matrix form as an effec- 
tive eigenvalue equation, which needs to be solved self- 
consistently^ 5 - 

^ H rr i(j) ST , = £s<Psr, (Al) 

where the effective Hamiltonian matrix is given by 

H rr < = h TT i + Mrr' ^«rr" \4>sr"\ 2 - ^rr'^sr'^r^ ■ 

(A2) 

The indices r, r' and r" indicate positions on the lattice, 
h rr r is the single-particle part of the Hamiltonian, and 
v rr r is the interaction energy of a particle at site r and a 
particle at site r'. cf> sr are the eigenvectors in lattice-site 
representation, each indexed by label s and with s s as 
its eigenvalue. The self-consistent solution yields the HF 
ground state energy, given by 

s— 1 \ rr' / 

where the summation is over the n states with lowest 
eigenvalues e s . For the Hamiltonian, Eq. ([3]), considered 
in the present study the single-particle term is 

*~'=¥-Wi^r (A4) 

with an implicit dependence, given by Eq. of t rT > 

on the configuration of the core spins {SV}- The sec- 
ond term in Eq. (| A4|) reflects the interaction between 



12 



the conduction electrons and a continuous neutralizing 
positive background of density x — n/A, where A is the 
system area, via the Coulomb potential v rr > = U/ |r — r'|. 
Noting that the eigenvectors are normalized to unity, 



J2r" \4>sr" | = lj and that the r = r" and r = r' terms in 
Eq. (|A2p are equal and opposite, we are led to analyze 
the following HF Hamiltonian 



s=l 



dR 
7 ~^R\ 



(A5) 



We are interested in cases where the core spins config- 
uration is periodic, such that the system can be divided 
into N unit cells, each containing an identical configura- 
tion of spins on M sites. Let the super-lattice vectors {1} 
identify the location of the unit cells. A position r on the 
lattice can then be written as r (I, i) = 1 + r,, where r$ 
is the position within the unit cell 1, containing r. The 
spin periodicity implies that if rr < between sites r = 1 + r 
and r' = 1' + rj depends only on i, j, and the super- 
lattice vector 1" = 1' — 1 connecting the two unit cells, 



i.e., H rr ' = Hij(\"). As a consequence of Bloch's theo- 
rem this means that the energy eigenf unctions, expressed 
in the (1, i) representation, take the form 0f,i(k)e ik l /\/]V, 
with eigenenergies £&(k), where k is defined within the 
first Brillouin zone of the reciprocal super-lattice. The 
"band" index b runs from 1 to M, and <fibi(k) is nor- 
malized to unity within a single unit cell. Written in 
the (k, i) basis, the Hamiltonian becomes block diago- 
nal, where the matrix elements of the block connecting 
states with the same k are given by 



J 



Hij (k) 



N ^ 

b,k' 



9 [/x - e b (k') 



6 ij Yl v " 



b w (kOr-%'(k-k')<&(k')0 w (k') 



(A6) 



Here 



Mk) = £Ml)e- 4k -\ 



(A7) 



is the Fourier transform of the hopping amplitudes tij (1) 
between sites i and j in unit cells separated by a super- 
lattice vector 1. We also introduced 



,.H 



(0) 



dR 



A u .c. J | r i R-l 



(A8) 
(A9) 



where A„. c . = A/N is the area of a unit cell. The chemi- 
cal potential, /i, is defined by n = J2b k £ b(k)], with 
0(x) denoting the step function. 

The HF ground-state energy of the conduction elec- 
trons is 



E, 



-E 

2 ^ 

b,k 



e b (k) + i^^(k)^(k)^(k) 



. 1 f dR 
8[/i-£ 6 (k)j- -Unx J — 



(A10) 



where in the last term we have taken the limit N — > 
oo. This diverging contribution is canceled by the self- 
interaction of the positive background, evaluated in the 
same limit 



dR dR' 1 
\R—W\ = 2 UnX 



w (A11) 



Consequently, the total energy (not including the contri- 
bution of the anti-ferromagnetic interaction between the 
core spins) is given by the sum over b and k in Eq. (|A10p . 

In order to evaluate the matrix elements of the HF 
Hamiltonian, we need a method to calculate the infinite 
super-lattice sums in Eqs. (|A7[) - (|A9[) . The first of these 
is trivial since hopping is allowed only between nearest- 
neighbor sites within a unit cell or between adjacent ones. 
In two dimensions this leaves at most five terms to the 
sum. On the other hand, the Coulomb interaction is 
long-ranged, and an infinite number of terms needs to be 
included in Eqs. (|A8]l and (fAljj) . 



The Hartree term 



The Hartree interaction matrix (|A9|1 includes two di- 
verging contributions, one coming from the interaction 
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with the average electronic density and the other from 
the interaction with the positive uniform background. 
The two contributions cancel each other. In order to 
demonstrate this and extract the remaining finite piece 
we employ Ewald summation (see Appendix [B]) . The 
main identity of this method, directly applicable to the 
evaluation of the first term in Eq. (|A9[) . is 



En 



1 



2tt 



- 1 , 



E 



1^ 6rfc V2(V 



ErfT^ erfc(G|1+r|) - ( A12 ) 



true for the Fock term since it satisfies 



TjFock 



(k) 



-E 

b,k' 



9(/i - e b )Vij (k - k') (j)* bi (j)bj, 



.[/^ec/i-^i— ^lifoj- (Ai5) 



Moreover, Eq. (|A15|) implies that in the case of FM 
droplets the calculation of the Fock term involves only a 
finite sum (over the M states within each droplet). This 
is a direct consequence of the vanishing overlap between 
electronic states in different droplets. 



As before, 1 are the super-lattice vectors and A u c , is the 
unit cell area. Here, g are the reciprocal super-lattice vec- 
tors, and G is an arbitrary constant, chosen to minimize 
the number of relevant terms in both sums controlled by 
the complementary error function erfc (x) . Note that the 
divergence which stems from summing over large 1 vec- 
tors in the left hand side of Eq. (|A12j) is encoded in the 
g = term on the right hand side. This divergence is 
canceled by the integral over the whole system in Eq. 
(IA9[) . This can be readily seen by using Eq. (|A12[) with 
G — > oo, to write it as 



1 fdr = 1 r dr sr_j_ 
l«.c. J |r| A u . c . J u , c , r V |! + 



r| A. 



2ir -^-^ S Si q 



Consequently we find for the Hartree matrix 



_ e*^ f\ S \\ 20F 

Au. c .2-i |g| 



2GJ Au.c.G 
X^erfcCGIl + ^l) 

1^0 

, crfc(G|r ri ,|) , 2G 
j j o U '—j=, (A14) 



where r,-,y = r, - r,/ 



2. The droplets Fock term 

When the core spins are arranged in FM droplets sep- 
arated by an AFM ordered background, the conduction 
electrons cannot hop from one unit cell to the other, i.e., 
tij (1) = tij5\ Q. Consequently Uj (k) is independent of 
k, see Eq. (|A7[) . Under such a condition it is easy 
to verify that the HF eigenfunctions and eigenenergies 
are k-independent as well. To prove this assertion, let 
us assume that it is true and show that it leads to a 
k-independent HF Hamiltonian, hence closing the argu- 
ment self-consistently. Since the Hartree term in the HF 
Hamiltonian, Eq. (|A6[) . depends on k only through the 
HF eigenfunctions it obviously fulfills the requirement. 
To complete the demonstration we note that the same is 



3. The stripes Fock term 

When the core spins are arranged in a striped config- 
uration, hopping is allowed between unit cells along the 
direction of the stripes. In other words, if we decompose 
the super-lattice vectors as 1 = n a a + n b h, where a and 
b are primitive vectors along and off the stripe direction, 
respectively, then tij (n a ,n b ) — Uj (n a ) S nbi Q. As a re- 
sult (k), depends only on the k component along the 
stripes, i.e., tij (k) = t^ (k a ). It follows then, using the 
same reasoning presented above for the droplet case, that 
the HF eigenfunctions and eigenenergies depend only on 
k a , and the Fock term takes the form 

U 



(A13) Hf^ k (k a ) = -—J2e^-e b (KW b dk' a )^(K) 



b,k' 



E 



1 - 5ij8n a S) in a (k a -k'„)a 



Tij + n a a\ 



(A16) 



where N a is the number of unit cells along the stripe, 
a = |a|, and r,j is the vector connecting sites i and j 
within a unit cell. 

In contrast to the Hartree term where the interac- 
tion decays slowly, the exponential factor in the Fock 
exchange, Eq. (|A16|) . ensures that the series converges 
relatively fast. Hence, the infinite sum is well approxi- 
mated by assuming a long but finite stripe. In our calcu- 
lation, we used N a = 100 — 200, and verified that larger 
values change the ground state energy by an insignificant 
amount. Note that the logarithmic divergence of the n a 
sum in the case k' a — k a , is integrable, and vanishes upon 
the summation over k' . 



APPENDIX B: EWALD SUMMATION IN 2D 



The development (based on Ref. 1561 ) of Ewald's sum- 
mation method begins with defining the function 



2 ^- e -|l+r|V 



(Bl) 



where the vectors {1} correspond to the N points of a two- 
dimensional lattice of area A. F is a periodic function of 



14 



r, with the periodicity of the lattice. Therefore, it can be 
expanded into the following Fourier series 



(B2) 



where {g} are the reciprocal lattice vectors, and 
2 



F S (P) 



|l+r|V e -igT 



2 N 

-|g| 2 /4p 2 



A-u.c.P 



(B3) 



Here, A u . c , = A/N is the area of a unit cell. Using Eqs. 
(|B2 ]) -(|B3 ]l and the identity 



1 



4= I'dpc 1 ■ 



(B4) 



we obtain 



- V- / - . e -|l+r|V 



«-c. g Jo P 



a. 



-|g| 2 /4/3 2 +ig-r 



X -.'e- 11 ^ 12 " 2 , (Br,) 

'7T 



1 



where the integral was split into two at an arbitrary pos- 
itive value G. Finally, calculating the integrals leads to 
Eq. (|A12|) . where the divergent piece of the original sum 
is given by the g = term (and when r = also the 1 = 
term) in the new representation. The remaining part of 
the infinite sums over g and 1 is rapidly converging at a 
rate which is optimized by an appropriate choice of G. 
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